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Abstract: We show that weakly guiding nonlinear waveguides support 
stable propagation of rotating spatial solitons (azimuthons). We investigate 
the role of waveguide symmetry on the soliton rotation. We find that 
azimuthons in circular waveguides always rotate rigidly during propagation 
and the analytically predicted rotation frequency is in excellent agreement 
with numerical simulations. On the other hand, azimuthons in square 
waveguides may experience spatial deformation during propagation. 
Moreover, we show that there is a critical value for the modulation depth 
of azimuthons above which solitons just wobble back and forth, and below 
which they rotate continuously. We explain these dynamics using the con- 
cept of energy difference between different orientations of the azimuthon. 
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1. Introduction 

Spatial solitons are nonlinear localized states that keep their form during propagation due to 
the balance between diffraction and self-induced nonlinear potential [1]. Recently, there has 
been a lot of interest in a generalized type of spatial solitons, the so-called azimuthons. These 
are azimuthally modulated beams, that exhibit steady angular rotation upon propagation [2]. 
They can be considered as azimuthally perturbed optical vortices, i.e. beams with singular 
phase structure [2-4] . Theoretical studies demonstrated both, stable and unstable propagation 
of azimuthons [5-7], and the first experimental observation of optical azimuthons was recently 
achieved in rubidium vapors [8]. 

It appears that higher order solitonic structures and optical vortices are generally unstable in 
typical nonlinear media with local (Kerr-like) response [9, 10]. On the other hand, it has been 
shown that a spatially nonlocal nonlinear response provides stabilization of various complex 
solitonic structures including vortices [11-14]. Consequently, azimuthons and their dynam- 
ics have been studied almost exclusively in the context of spatially nonlocal nonlinear me- 
dia [7, 15-17]. In spite of the fact that there are various physical settings exhibiting nonlocality 



such as nematic liquid crystals [18-20], Bose-Einstein condensates [4,21,22], plasmas [23], 
thermo-optical materials [24] etc., experimental realization of such systems is always quite in- 
volved. Moreover, from the theoretical point of view, nonlocal media are quite challenging for 
numerical modeling and analytical treatment. 

In this paper we propose a much simpler and experimentally accessible optical system to 
study the propagation of azimuthons: a weakly nonlinear optical multi-mode waveguide. Here, 
weakly nonlinear means that the nonlinear induced index change, which is proportional to the 
intensity of the optical beam, is small compared to the index profile (or trapping potential) of 
the waveguide. Following [25], we can expect that weakly nonlinear azimuthons are stable in 
multi-mode waveguides. 

The paper is organized as follows: In Sec. 2, we briefly introduce the general model equation 
for beam propagation in weakly nonlinear waveguides, and then we discuss in detail the prop- 
erties of dipole azimuthons in circular and square waveguides in Sec. 3, respectively. In Sec. 4, 
higher order azimuthons are investigated, and in Sec. 5 we conclude. 

2. Mathematical modeling 

We consider the evolution of a continuous wave (CW) optical beam with amplitude £ , 7] , Q, 
where (§,T]) and £ denote the transverse and longitudinal coordinates, respectively. Then the 
propagation of this beam in a weakly-guiding waveguide with Kerr nonlinearity in the scalar, 
slowly varying envelope approximation is described by the following equation [26] : 

where Icq = 27r^/Ao refers to the carrier central wave number in the medium; is the Kerr 
nonlinear coefficient, n(§ , rj) the linear refractive index distribution, rib the background index, 
and lim^ ^^^(£,77) = n^. We consider a weakly-guiding waveguide, so both the linear and 
nonlinear induced index change \n — rib\ and ri2 1 £\ 2 are small compared to the mean index rib, 
and at the same time ri2\£\ 2 <C \n — rib\, to guarantee a weak nonlinearity. 

From the mathematical point of view, Eq. (1) is a (2+l)-dimensional nonlinear Schrodinger 
(NLS) equation with linear potential representing the waveguide profile. For technical conve- 
nience, we rescale Eq. (1) to dimensionless quantities with x = y = Tf/ro, z = £/(2&o?o), 
and G = sgn(ft2), where ro represents the transverse spatial extent of the waveguide. Then, the 
two-dimensional (2D) NLS equation for the scaled wave function y/ = koroy/2\ri2\/ni ) £ reads 

with an "attractive" bounded potential V = 2k^r\(n — rib) /rib, given by the spatial profile of the 
refractive index of the waveguide. Here we will consider propagation in step index waveguides 
with parameter V(x,y) given as follows 



V(x,y) 



Vq where 
elsewhere 



\Jx 2 + y 2 < 1 for circular waveguide 
|jc| < 1 & \y\ < 1 for square waveguide 



Vb is the height of the potential and determines the number of linear waveguide modes. 

Considering the specific case of silica-made waveguides, typical values for the parame- 
ters are rib = 1-4, \n — rib\ < 9 x 10 -3 , ^ = 3x 10 -16 cm 2 /W at a vacuum wavelength of 
Aq = 790 nm, and those values do not vary much when we choose any Aq in the range from 



visible to near-infrared. The laser-induced-damage-threshold (LIDT) for synthetic fused sil- 
ica is about 10 J/cm 2 for nanosecond pulses (according to Eq. (1) and Table 2 in Ref. [27]) 
corresponding to ~ 10 10 W/cm 2 , up to which the model should be valid. For shorter, femtosec- 
ond pulses this damage threshold is 1000 times higher, but since we consider "CW" beams 
the lower value for nanosecond pulses is relevant. It is justified to neglect temporal effects 
on beam propagation, because the dispersion length is already of the order of kilometers for 
pulse durations of a few tens of picoseconds, much longer than propagation distances consid- 
ered in this work. Throughout this paper, we choose Vb = 1000, and find that ro « 25.0 jim 

through ro = \J~\n^V~/ [2&q(^ — n^)] |, which nicely corresponds to the radius of a standard cir- 
cular multi-mode fiber [28] . In addition, in order to guarantee intensities smaller than the LIDT 
the condition | \\f\ 2 < 1 /3 should be fulfilled. 

Note that in Eq. (2) a is the sign of nonlinearity: a = 1 (a = — 1) represents focusing (defo- 
cusing) nonlinearity. The main difference between weakly nonlinear waveguides with focusing 
and defocusing nonlinearity is that defocusing nonlinearity supports higher amplitudes of the 
wave function, and in general leads to more stable and robust configurations. In this paper, 
however, we consider the experimentally relevant case of a focusing nonlinearity. 

3. Rotating localized dipoles 

3.1. Dipole azimuthons in circular waveguides 

Azimuthons are a straightforward generalization of the usual ansatz for stationary solutions 
[2] . They represent spatially rotating structures and hence involve an additional parameter, the 
rotation frequency CO (see also [29]), so we seek approximate solutions of the form 

\lf(r,<l>,z)=U(r,<l>-m)exp(itcz), (3) 

where r = ^Jx 2 + y 2 and (f) the azimuthal angle in the transverse plane (x,y), U is the stationary 
profile, CO the rotation frequency, and K the propagation constant. For (0 = 0, azimuthons be- 
come ordinary (non-rotating) solitons. The simplest family of azimuthons is the one connecting 
the dipole soliton with the single charged vortex soliton [15]. A single charged vortex consists 
of two equal-amplitude dipole-shaped structures representing real and imaginary part of U. If 
these two components differ in amplitude, the resulting structure forms a "rotating dipole" az- 
imuthon. If one of the components is zero we deal with the dipole soliton, which consists of 
two out-of-phase humps and does not rotate for symmetry reasons. In a first attempt, let us as- 
sume we know the radial shape of the linear vortex mode F(r) which is normalized according 
to 71 J r\F(r)\ 2 dr = 1. Then, using separation of variables, we consider the simplest so-called 
"rotating dipole" azimuthon with ansatz [16] 

£/(r, (j) - coz) = AF(r) [cos(0 - coz) + /£sin(0 - coz)} , (4) 

where A is an amplitude factor, and 1 — B the azimuthal modulation depth of the resulting ring- 
like structure. Because we are operating in the weakly nonlinear regime, using linear waveguide 
modes as initial conditions for nonlinear (soliton) solutions is a quite good approximation. 

After plugging Eq. (3) into Eq. (2), multiplying by U* and dU*/d(j> respectively, and inte- 
grating over the transverse coordinates we end up with [16] 

- kP+ coL z + I + N = 0, (5a) 



-kL z + cdP , + I , + N' = 0. 



(5b) 



This system relates the propagation constant K and the rotation frequency CO of the azimuthons 
to integrals over their stationary amplitude profiles: 



P = J J r\U(r)\ 2 drd$, (6a) 

L z = -ijj r^V(r) dr (6b) 

/ = jj rU* (r) A ± U (r) dr d0 , (6c) 

= jjr[o\U(r)\ 2 + V]\U(r)\ 2 drd$, (6d) 



'-//- 



I' = ijj r^^A ± U(r) dr d(j>, 
N' = iJJ r[a\U(r)\ 2 + V] ^^-U{r) drdf 



drd</>, (6e) 

(6f) 



(6g) 



The first two quantities (P and L z ) have straightforward physical meanings, namely power 
and angular momentum. The integrals / and /' are related to the diffraction mechanism of the 
system, whereas N and N' account for waveguide and nonlinearity. Thus we can formally solve 
for the rotation frequency with these quantities and obtain 

P(I>+N')-L Z {I + N) 

(0 = — • () 

After inserting Eq. (4) into Eq. (7), it turns out that only the nonlinear term contributes to the 
rotation frequency CO (see also Eq. (10) in Ref. [30]). In order to give an estimate for the rotation 
frequency, we use the linear stationary modal profile F(r) of a circular waveguide expressed in 
terms of Bessel functions of first kind J\ and modified Bessel function of second kind K\ : 



{Ji(VVo-Kr) for0<r<l 
Ji(y/Vo=K) „ , r- v f ^ , i (8) 
Mft l(v ^ r) 

where C is a normalization factor such that 71 J \F(r)\ 2 rdr = 1. Thus, the analytically predicted 
rotation frequency is 

CO =™ j r\F(r)\ 4 dr-A 2 B. (9) 

From the above equation, one can see that sense of rotation is opposite for focusing and defo- 
cusing nonlinearities. Moreover, as CO = 2&o?o -2n/£, one can find from Eq. (9) an expression 
for the physical distance £ over which the azimuthon in a circular waveguide performs one full 
rotation: 

8 *° r ° 2 (10) 



A 2 B J r\F(r) 



For example, assuming propagation in a silica fiber with the parameters discussed earlier (see 
Sec. 2) and taking A = 0.4, B = 0.5, we obtain £ « 2.4 m. 



Figure 1 shows the dependence of the azimuthon rotation frequency as a function of its 
amplitude A (left panel) and the modulation parameter B (right panel). Blue solid lines represent 
the analytical predictions and, in excellent agreement, red dots are obtained from numerical 
simulations. An exemplary propagation of the dipole azimuthon in circular waveguide rotating 
with CO w 0.037 is illustrated in Fig. 2. The top iso-surface plot displays 3D rigid and continuous 
rotation of the azimuthon during propagation. The row below depicts the transverse intensity 
(left) and phase (right) distribution of the dipole azimuthon during propagation after rotating 
by x/2. 




Fig. 1. Azimuthon rotation frequency CO versus amplitude factor A (left panel) and ampli- 
tude ratio B (right panel). Blue solid lines show analytical predictions from Eq. (9). Red 
dots denote results obtained from numerical simulations of Eq. (2). The values of A and B 
are shown next to the lines. 
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Fig. 2. The propagation of a dipole azimuthon with A = 0.4, B = 0.5 in the circular waveg- 
uide. The iso-surface plot on the top clearly displays the spiraling of the azimuthon during 
propagation. Figures in the row below depict the azimuthon' s transverse intensity distribu- 
tion at input, and after rotating n/2, respectively. The white line indicates the waveguide 
boundaries. 



3.2. Azimuthon-like dipoles in square waveguides 

Because of the lack of circular symmetry, azimuthons in the strict sense of Eq. (4) (preservation 
of shape during rotation with constant frequency) cannot exist in a square waveguide. However, 
as we will show below, the azimuthon-like behavior is still possible even though the beam 
propagation is accompanied by variation of the beam transverse intensity distribution. To set 
the initial condition, we use two linear degenerated orthogonal dipole modes D\, D2 (as in the 
case of the circular waveguide), which are normalized according to JJ \D\^\ 2 dx dy = I, and 
superpose them as before to form the azimuthon-like object 

U{xmz = 0) = A (Di + iBD 2 ) • (11) 

The field Eq. (11) is then used as an initial condition to the nonlinear Schrodinger equation 
Eq. (2). However, in contrast to the circular waveguide, the orientation of the dipoles D\ and 
Z>2 is important in nonlinear regime. If the two orthogonal dipoles D\ and D2 are oriented along 
the diagonals of the waveguide cross-section (see first subplot of Fig. 3), for a given amplitude 
A rotation occurs only if the modulation parameter B exceeds a critical value B = B cr . Moreover, 
the rotation is no longer constant as in the case of cylindrical waveguide but fluctuates and hence 
in what follows we use its average value (termed "average frequency") (b = 1/LJq CO(z) dz with 
L being the propagation distance corresponding to one full 2n rotation. This threshold value B CY 
decreases if we choose different initial dipole orientations 1 , and appears to be zero (within our 
numerical accuracy) for parallel orientation with respect to the waveguide boundaries. 



INTENSITY z=0.00 PHASE INTENSITY z=30.00 PHASE 




PROPAGATION DISTANCE 



Fig. 3. The propagation of the azimuthon-like dipole, that rotates in the square waveguide. 
Top row: intensity (left) and phase (right) distribution corresponding to the input and soliton 
rotation by tf/4, respectively. The white line indicates the waveguide boundaries. The iso- 
surface plot at the bottom displays the rotation and deformation during propagation. The 
initial amplitude and modulation parameters are A = 0.4, B = 0. 5. 

Let us focus on the case where the initial dipole-like field structure is oriented along the 
diagonal of the waveguide cross-section, and from now on the notation D\ and D2 stands for this 
orientation. Figure 3 shows an exemplary propagation of the resulting azimuthon-like solution 



Tn the linear system, any dipole orientation is possible, and can be constructed from orthogonal basis (Z>i,Z>2) by 
superposition. 



with B > B cr , and an average rotation frequency of G) « 0.0262. The top row depicts intensity 
and phase distribution at different propagation distance while the bottom plots illustrate the 
full 3D evolution of the "soliton". The rotation (counter-rotating w.r.t. phase) accompanied 
by beam deformation is clearly visible. For an amplitude ratio B less than the critical value 
(B < B cr ) the azimuthon no longer rotates but swings back and forth and hence its average 
rotation frequency is zero. The propagation of this wobbling dipole with an amplitude ratio B 
smaller than the critical value is displayed in Fig. 4. The right panels show the maximum angle 
which the dipole attains during propagation. The 3D surface plot in the bottom row clearly 
illustrates the swinging movement of the "soliton" upon propagation. 
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Fig. 4. The propagation of a dipole azimuthon that wobbles in the square waveguide. The 
first row depicts the dipole at input and maximum displacement, respectively. The white 
line indicates the waveguide boundaries. The iso- surface plot below displays the twist and 
deformation during propagation. The chosen amplitude factor and ratio are A = 0.4, B = 
0.2. 

Due to the azimuthon profile deformation, it is not possible to find an analytical expression 
of d) as a function of A and B as in the case of circular waveguide (i.e. Eq. (9)). Therefore, we 
need to resort to numerical simulations. In Fig. 5 we show the numerically determined relation 
between cb and the azimuthon parameters A (left panel) and B (right panel). The threshold-like 
behavior is evident in the right graph of this figure: the region B < B cr , in which the soliton 
wobbles, is depicted by a green line. 

Let us have a closer look at the boundary between domains of azimuthon rotation (B > B cr ) 
and wobbling (B < B cr ). We find numerically that the critical value B CY , which separates those 
two domains, depends very weakly on the azimuthon amplitude A: we are in weakly nonlinear 
limit, and may use linear waveguide modes to elucidate nonlinear propagation properties of the 
azimuthons. The appearance of a threshold value B cr , above which the azimuthon rotates, can 
be explained by considering the Hamiltonian associated with propagation equation Eq. (2), 

II (\^¥\ 2 -hw\ 4 -V\w\ 2 ) d*dy, (12) 



which is a conserved quantity upon propagation. Depending on their orientation, stationary 
dipoles (B = 0) have a slightly different value of Jt? in nonlinear regime. If we take our diagonal 



0.06 




Fig. 5. Averaged rotation frequency co versus amplitude factor A (left panel) and amplitude 
ratio B (right panel) in a square waveguide. Red dots are numerical results obtained from 
Eq. (2), blue curves (B < 0.4) and blue lines {B > 0.4) are the fitting results to the numerical 
simulations. The green line represents the superposed dipoles twist during propagation. The 
values of A and B are shown next to the curves. The black square represents the analytical 
estimate for B CY . 



(w.r.t. waveguide cross-section) dipoles D\ and D 2 from above, we can construct a parallel 
dipole D p in the following way: 



jj \D l ^D 2 \ 2 dxdy 

The intensity distribution of this dipole D p is very close to that shown in the snapshot of the 
rotating azimuthon shown in the right panel in Fig. 3. One can show using the definition of 
the Hamiltonian Eq. (12) that J^(D p ) > Jt?{D x ) = Jf?(D 2 ). In fact, it turns out that those 
two dipole orientations (parallel and diagonal) correspond to the extremal values of the Hamil- 
tonian. Now, if the azimuthon rotates it has to pass through all possible orientations, i.e., its 
value of has to be larger than Jt?(D p ). We can use this reasoning to estimate the value of 
B cr : The left panel of Fig. 6 depicts the dependence of the Hamiltonian of superposed diag- 
onal dipoles as a function of B (red line) for constant power P. The blue line represents the 
value of the Hamiltonian of the parallel dipole ~ D p at the same power level. As the red curve 
monotonically increases with modulation parameter B, the intersection of the two curves gives 
an estimate of the critical value of the azimuthon modulation B cr , which is shown in the right 
panel in Fig. 5 by a black square. In other words, in order to rotate from the diagonal to vertical 
position the azimuthon must overcome some kind of energy barrier represented by the differ- 
ence of the values of the Hamiltonian in those two basic states. As Fig. 6 shows this can be 
achieved by increasing the value of B in the diagonal state to B\ = 0.342, which is very close 
to the value found numerically (see right panel in Fig. 5). Similar phenomena have been identi- 
fied, for instance, in discrete nonlinear systems where the discreteness introduces the so called 
Peierls-Nabarro potential which has to be overcome by a soliton to become mobile [31]. 

It should be emphasized that the crucial difference between azimuthon dynamics in circular 
and square waveguide originates from the fact that angular momentum of the beam is conserved 
in the circular waveguide, whereas it changes considerably in the square waveguide. In the right 
panel of Fig. 6 we display the angular momentum of the beam in the square waveguide for 
four different values of modulation parameter B. If B is above the critical value B cr , the angular 
momentum does not change its sign as the blue curve shows. If B < B cv , the angular momentum 
changes periodically its sign in propagation (red curve), which indicates the swinging motion 
of the soliton. For values B close to the critical value the variation of the angular momentum 
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Fig. 6. Left panel: Hamiltonian of A{D\ + iBD2) (red curve) as a function of B and Hamil- 
tonian of D3 (blue line) (for the same power P). The intersection of the two curves defines 
B CY (depicted by the black square). Right panel: The angular momenta of the rotating super- 
posed dipoles (corresponding to Fig. 3, blue curve), and the twisting ones (corresponding to 
Fig. 4, red curve). The other curves depict the cases of B slightly above (black) and bellow 
(green) B CY . All curves are obtained for A = 0.4. 



are largest as shown by the green and black curves. At the same time, the appearance of broad 
plateaus indicates slow propagation dynamics. 

It is worth mentioning that we also investigated other waveguide symmetries, such as hexag- 
onal waveguides, and found the behavior of the dipole azimuthons to be qualitatively similar to 
that of the square waveguide, i.e., they either rotate or wobble, depending on modulation depth 
and initial orientation. However, the difference of the values of the Hamiltonian for the two 
extremal dipoles is now smaller, resulting in a lower threshold value B cr . We believe that this 
threshold behavior in the propagation dynamics of azimuthons is generic for any non-circular- 
symmetric waveguide structure which supports degenerated dipole modes. 

4. Rotating higher order localized modes 

In analogy to dipole-azimuthons discussed so far, it is possible to consider dynamics of 
higher order azimuthons constructed by using pairs of degenerated higher order waveguide 
modes. In particular, in a circular waveguide degenerate higher order modes have the form of 
quadrupoles, hexapoles, octapoles, decapoles, etc. and other modes with even number of lobes 
(optical necklaces). In a square waveguide, one can identify pairs of degenerate modes in a 
form of hexapoles, octapoles, decapoles, dodecapoles, etc. (optical matrices). Interestingly, the 
quadrupoles found in square waveguides are not degenerate and hence the propagation dynam- 
ics of their corresponding superposed states in weakly nonlinear regime is dominated by linear 
mode beating. 

Fig. 7 displays an example of hexapole azimuthon in a circular waveguide. The left panel 
represents the initial intensity and phase distribution of the azimuthon. The 3D surface plot 
in the right panel depicts stable and smooth rotation of the hexapole azimuthon. As far as 
the higher order azimuthons of the square waveguide are concerned, they share the dynamical 
properties of dipole azimuthons discussed before. As an example, Fig. 8 shows the rotating (top 
row) and swinging (middle row) hexapoles. Due to the beam deformation, the hexapole with 
B = 0.5 transforms into a (2 x 3) matrix while rotating. The left two panels in the third row of 
Fig. 8 display the rotating dodecapole azimuthon and its deformation into a (3 x 4) solitonic 
matrix; the right panels show rotating icosapole azimuthon and its deformation into a (4 x 5) 
solitonic matrix. 




Fig. 7. The propagation of a hexapole azimuthon with A = 0.4, B = 0.5 in the circular 
waveguide. Left panel shows the original hexapole azimuthon and the corresponding phase; 
right panel shows the iso- surface plot of the propagation. 
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Fig. 8. The propagation of superposed hexapoles in the square waveguide. The first row 
shows a rotating azimuthon with A = 0.4, B = 0.5 > B CY (Media 1), and the second row 
shows a twisting azimuthon with A = 0.4, B = 0.2 < B CY (Media 2). The superposed do- 
decapoles (left two panels), and superposed icosapoles (right two panels) as well as their 
deformations are shown in the third row; in both cases we choose A = 0.4 and B = 0.5. 



5. Conclusion 



In conclusion, we demonstrated numerically stable propagation of azimuthons, i.e. localized 
rotating nonlinear waves in weakly nonlinear waveguides. Depending on the waveguide pro- 
file, different nonlinearity induced propagation dynamics can be observed. We showed that az- 
imuthons in circular waveguides rotate continuously. The analytically predicted dependence of 
rotation frequency ft) as a function of soliton parameters was found to be in excellent agreement 
with numerical simulations. Further on, we discussed propagation of azimuthon-like structures 
in square waveguides. We showed that their dynamics critically depends on the initial condi- 
tions. In particular, we found a threshold-like behavior in the propagation dynamics, separating 
rotating and wobbling solutions. We showed that this effect is associated with different values 
of the Hamiltonian for different azimuthon orientations. As our analysis relies on physical pa- 
rameters of actual multi-mode waveguides, our findings may open a relatively easy route to 
experimental observations of stable azimuthons. 
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